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Abstract 

Nanomanipulation of individual RNA molecules, using laser optical tweezers, has made it pos- 
sible to infer the major features of their energy landscape. Time dependent mechanical unfolding 
trajectories, measured at a constant stretching force (fs), of simple RNA structures (hairpins and 
three helix junctions) sandwiched between RNA/DNA hybrid handles show that they unfold in a 
reversible all-or-none manner. In order to provide a molecular interpretation of the experiments 
we use a general coarse-grained off-lattice Go-like model, in which each nucleotide is represented 
using three interaction sites. Using coarse-grained model we have explored forced-unfolding of RNA 
hairpin as a function of fs and the loading rate (rj). The simulations and theoretical analysis have 
been done without and with the handles that are explicitly modeled by semiflexible polymer chains. 
The mechanisms and time scales for denaturation by temperature jump and mechanical unfolding 
are vastly different. The directed perturbation of the native state by fs results in a sequential un- 
folding of the hairpin starting from their ends whereas thermal denaturation occurs stochastically. 
From the dependence of the unfolding rates on rj and fs we show that the position of the unfolding 
transition state (TS) is not a constant but moves dramatically as either rj or fs is changed. The 
TS movements are interpreted by adopting the Hammond postulate for forced-unfolding. Forced- 
unfolding simulations of RNA, with handles attached to the two ends, show that the value of the 
unfolding force increases (especially at high pulling speeds) as the length of the handles increases. 
The pathways for refolding of RNA from stretched initial conformation, upon quenching fs to the 
quench force /q, are highly heterogeneous. The refolding times, upon force quench, are at least an 
order of magnitude greater than those obtained by temperature quench. The long /g-dependent 
refolding times starting from fully stretched states are analyzed using a model that accounts for 
the microscopic steps in the rate limiting step which involves the trans to gauche transitions of 
the dihedral angles in the GAAA tetraloop. The simulations with explicit molecular model for the 
handles show that the dynamics of force-quench refolding is strongly dependent on the interplay 
of their contour length and the persistence length, and the RNA persistence length. Using the 
generality of our results we also make a number of precise experimentally testable predictions. 
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2 



INTRODUCTION 



RNA molecules adopt precisely denned three dimensional structures in order to perform 
specific functions 0|. To reveal the folding pathways navigated by RNA en route to their 
native conformations requires exhaustive exploration of the complex underlying energy land- 
scape over a wide range of external conditions. In recent years, mechanical force has been 
used to probe the unfolding of a number of RNA molecules 0, HI, 0] • Force is a novel way 
of probing regions of the energy landscape that cannot be accessed by conventional meth- 
ods (temperature changes or variations in counterion concentrations). In addition, response 
of RNA to force is relevant in a number of cellular processes such as mRNA translocation 
through the ribosome and the activity of RNA-dependent RNA polymerases. Indeed, many 
dynamical processes are controlled by deformation of biomolecules by mechanical force. 

By exploiting the ability of single molecule laser optical tweezer (LOT) to control the 
magnitude of the applied force Bustamante and coworkers have generated mechanical un- 
folding trajectories for RNA hairpins and Tetrahymena thermophila ribozyme 0|. The 
unfolding of the ribozyme shows multiple routes with great heterogeneity in the unfolding 
pathways 0. In their first study [H[ they showed that simple RNA constructs (P5ab RNA 
hairpins or a three helix junction) unfold reversibly at equilibrium. From the time traces 
of the end-to-end distance (R) of P5ab, for a number of force values, Liphardt et. al. (H 
showed that the hairpins unfold in a two-state manner. The histograms of time dependent 
R (and assuming ergodicity) were used to calculate the free energy difference between the 
folded and unfolded states. Unfolding kinetics as a function of the stretching force fs was 
used to identify the position of the transition states 0, 0, @j- These experiments and sub- 
sequent studies have established force as a viable way of quantitatively probing the RNA 
energy landscape with R serving as a suitable reaction coordinate. 

The experiments by Bustamante and coworkers have led to a number of theoretical and 
computational studies using a variety of different methods 0, EH U, El, EH 3|- These 



studies have provided additional insights into the mechanical unfolding of RNA hairpins and 
ribozymes. In this paper we build on our previous work [l4[ and new theoretical analysis 
to address a number of questions that pertain to mechanical unfolding of RNA hairpins. In 
addition, we also provide the first report on force-quench refolding of RNA hairpins. The 
present paper addresses the following major questions: 

1) Are there differences in the mechanism of thermal and mechanical unfolding? We expect 

these two processes to proceed by different pathways because denaturation induced 
by temperature jump results in a stochastic perturbation of the native state while 
destabilization by force is due to a directed perturbation. The molecular model for 
P5GA gives a microscopic picture of these profound differences. 

2) For a given sequence does the position of the transition state move in response to changes 

in the loading rate (r/) or the stretching force /g? Based on the analysis of experi- 
ments over a narrow range of conditions (fixed temperature and loading rate) it has 
been suggested that the location of the sequence-dependent unfolding transition state 
(TS) for secondary structure is midway between the folded state while the TS for the 
ribozyme is close to the native conformation (Hf. Explicit simulations show that the 
TS moves dramatically especially at high values of fs and rj. As a consequence the 
dependence of the unfolding rate on fs deviates from the predictions of the Bell model 
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3) What is the origin of the dramatic differences in refolding by force-quench from stretched 

conformations and temperature-quench refolding? Experiments by Fernandez and Li 
on refolding initiated by force-quench on polyubiquitin construct suggest similar 
differences in the refolding time scales. For RNA hairpins we show that the incompat- 
ibility of the local loop structures in the stretched state and the folded conformations 
lead to extremely long refolding times upon force-quench. 

4) What are the effects of linker dynamics on forced-unfolding and force-quench refolding 

of RNA hairpins? The manipulation of RNA is done by attaching a handle or linker 
to the 3' and 5' ends of RNA. Linkers in the LOT experiments, that are done under 
near-equilibrium conditions, are RNA/DNA hybrid handles. These are appropriately 
modeled as worm-like chain (WLC). By adopting an explicit polymer model for semi- 
flexible chains we show that, under certain circumstance, non-equilibrium response 
of the handle (which is not relevant for the LOT experiments) can alter the forced- 
unfolding dynamics of RNA. We probe the effects of varying the linker characteristics 
on the dynamics of folding/ unfolding of RNA for the model of the handle- RNA-handle 
construct. In certain range of Tf non-equilibrium effects on the dynamics of linkers can 
affect the force-extension profiles. 

METHODS 

Hairpin sequence: To probe the dynamics of unfolding and force-quench refolding 
we have studied in detail the 22-nucleotide hairpin, P5GA whose solution structure has 
been determined by NMR (Protein Data Bank (PDB) id:leor). In many respects P5GA is 
similar to P5ab in the P5abc domain of group I intron 0] . Both these structures have GA 
mismatches and are characterized by the presence of the GAAA tetraloop. The sequence of 
P5GA is GGCGAAGUCGAAAGAUGGCGCC [TJ. 

RNA model: Because it is difficult to explore, using all atom models of biomolecules 
in explicit water, unfolding and refolding of RNA hairpins over a wide range of external 
conditions (temperature (T), stretching force (fs), and quench force (/q)), we have 
introduced a minimal off-lattice coarse-grained model 0] (Throughout the paper we use / 
for referring to mechanical force in general terms while fs and /q have specific meaning) . In 
this model each nucleotide is represented by three beads with interaction sites corresponding 
to the phosphate (P) group, the ribose (S) group, and the base (B) Thus, the RNA 

backbone is reduced to the polymeric structure [(P — S) n ] with the base that is covalently 
attached to the ribose center. In the minimal model the RNA molecule with N nucleotides 
corresponds to 3N interaction centers. The secondary structure and the lowest energy 
structure using the minimal model are shown in Fig^ 

Energy function: The total energy of a RNA conformation, that is specified by the 
coordinates of the 3N sites, is written as Vr = Vbl + Vba + Vdih + Vstack + Vnon + Velec- 
Harmonic potentials are used to enforce structural connectivity and backbone rigidity. The 
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connectivity between two beads (PjSj, SjPi+i and P^Si) is maintained using 

2N-2 N 

v BL = E M\nsp) i+1 - nsp),\ - (R° S p)i} 2 + E 2 kr ^ ~ fs > 1 ~ {R ° Bs)l}2 (1) 

1=1 8=1 

where fc r = 20 kcal/(mol ■ A 2 ), (R° SP )i and (R°Bs)i are the distances between covalently 
bonded beads in PDB structure. The notation (SP)i, denotes the i th backbone bead S or 
P. The angle 9 formed between three successive beads (Pi — Si — Pi + \ or Si-i — P, — S^) 
along sugar-phosphate backbone is subject to the bond-angle potential, 



v^ = E 2*»^-^ a (2) 



i=l 



where fcg = 20 kcaljimol ■ rad 2 ), and 6*° is the value in the PDB structure. 

Dihedral angle potential : The dihedral angle potential (Vdih) describes the ease of ro- 
tation around the angle formed between four successive beads along the sugar-phosphate 
backbone (^-ipSiP+i or PiSiPi+iSi+i). The z-th dihedral angle fa, which is the angle 
formed between the two planes defined by four successive beads % to % + 3, is defined by 
cosipi = (r*j + i 5 j x r i+ i t i +2 ) ■ (ri+2,t+i X ^+2,1+3)- In the coarse-grained model the right-handed 
chirality of RNA is realized by appropriate choices of the parameters in the dihedral po- 
tential. Based on the angles in the PDB structure (0°), one of the three types of dihedral 
potentials, trans (t, < <p° < 2tt/3), gauche(+) (g + , 2n/3 < <fi° < 47r/3), gauche(— ) (g~, 
4tt/3 < (f)° < 2ir), is assigned to each of the four successive beads along the backbone. The 
total dihedral potential of the hairpin is 

2N-4 

v DIH = E Mi + B l + c l 

i=X 

+ Alcos{^ - <j>° + + P 2 >s3(0, - # + + C* *n(& - 0° + #)] (3) 

where the parameters (in kcal/mol) defined for t, g + , and g~ are 

A u = 1.0, A 2i = -1.0, B u = B 2i = 1.6, Cu = 2.0, C 2i = -2.0 (77 = g+), 
A u = 1.0, A 2i = -1.0, B u = B 2l = 1.6, C u = 2.0, C 2i = 2.0 (77 = g~), 
A u = A 2i = 1.2, B u = B 2i = 1.2, C u = C 2i = 0.0 (77 = t). 

To account for the flexibility in the loop region we reduce the dihedral angle barrier by 

halving the parameter values in 19 < i < 24. 

Stacking interactions: Simple RNA secondary structures, such as hairpins, are lar gely 
stabilized by stacking interactions whose context dependent values are known 0, 0, 2(j . 



The folded P5GA RNA hairpin is stabilized by nine hydrogen bonds between the base pairs 
(see Fig^-B) including two GA mismatch pairs [T3| • The stacking interactions that stabilize 
a hairpin can be written as Vstack = Y^i=i x ^ ( n max = 8 in P5GA). We assume that the 
orientation-dependent Vi is 

Vi({(f>}, {tp}, {r}; T) — AGAT) x e - a st{sm 2 (0ii-</>?J+sm 2 (</> 2i -<^^ 

x e-Minj-r^+in+ij-i-r^) 2 } 
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where AG(T) = AH — TAS, the bond angles {0} are <f>u = ZSiBiBj, 2 « = ZBiBjSj, 
(p 3i = ZS i+ iB i+ iBj_i, </> 4 j = ZB i+1 Bj_iSj-i, the distance between two paired bases 
Tij = \Bi — Bj\, r i+ ij_i = \B i+ i — Bj-i\, and ipu and ip 2 i are the dihedral angles formed by 
the four beads BiSiSi + iB i+ i and Bj_iSj-iSjBj, respectively. The superscript o refers to 
angles and distances in the PDB structure. The values of a s t, f3 s t and 7^ are 1.0, 0.3A -2 
and 1.0 respectively. We take AH and AS from Turner's thermodynamic data set [lil lit . 
There are no estimates for GA related stacking interactions. Nucleotides G and A do not 
typically form a stable bond and hence GA pairing is considered a mismatch. We use the 
energy associated with GU for GA base pair. 

Nonbonded interactions: We use the Lennard- Jones interactions between non-bonded 
interaction centers to account for the hydrophobicity of the purine/pyrimidine groups. The 
total nonbonded potential is 

N-l N N 2N-1 2JV-4 2JV-1 

V NON = J2Y1 VBiBjir) + J2Y1 ' V B.(SP) m (r) V (SPU(SP) n (r) (5) 

i=l j=i+l i=l m=l m=l n=m+3 

where r = |r*j — fj\. The prime in the second term on the Eq.fjSJ denotes the condition 
m ^ 2i — 1. In our model, a native contact exists between two non-covalently bound beads 
provided they are within a cut-off distance r c (=7.0A). Two beads beyond r c are considered 
to be non-native. For a native contact, 

W0=^*[(^) U -2(^) 6 ] (6) 

where r°- is the distance between beads in PDB structure and Cf l lVj = 1.5 kcal/mol for all 
native contact pairs except for B w Bi 3 base pair associated with the formation of the hairpin 
loop, for which cj^ loBl3 = 3.0 kcal/mol. The additional stability for the base pair associated 
with loop formation is similar to the Turner's thermodynamic rule for the free energy gain 
in the tetraloop region. For beads beyond r c the interaction is 

Vto{r) = C R [fy" + fy*] (7) 

with a = 3.4A and Cr = 1 kcal/mol. The value of C^ Vj (= l.hkcal/mol) has been chosen 
so that the hairpin undergoes a first order transition from unfolded states. Our results are 
not sensitive to minor variations in CV 

Electrostatic interactions: The charges on the phosphate groups are efficiently screened 
by counterions so that in the folded state the destabilizing contribution of the electrostatic 
potential is relatively small. However, the nature of the RNA conformation (especially 
tertiary interactions) can be modulated by changing counterions. For simplicity, we as- 
sume that the electrostatic potential between the phosphate groups is pairwise additive 
Velec = J2i=~i 1 J2f=i+i V P t Pj( r )- For VPiPjir) we use Debye-Hiickel potential, which ac- 
counts for screening by condensed counterions and hydration effects, and is given by 

zp.zp.e 2 ,, 

V PiP] = ^3— e -rno 8 
3 4ne e r r 
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where Zp. = — 1 is the charge on the phosphate ion, e r = e/e and the Debye length 



1 d = y 8^**^1 witn ^ eiec = = 8 - 99 x 10 9 JmC~ 2 . To calculate the ionic strength 
/ = l/2'^2 i zfci, we use the value q = 200mM-NaCl from the header of PDB file 
[]~7l | . We use e r = 10 in the simulation j2jj. Because the Debye screening length ~ VT 
the strength of electrostatic interactions between the phosphate groups are tempera- 
ture dependent even when we ignore the variations of e with T. At room temperature 
(T ~ 300 K) the electrostatic repulsion between the phosphate groups at r ~5.8 A, 
which is the closest distance between phosphate groups, is Vp i p i+1 ~ 0.5 kcal/mol. Thus, 
Velec between phosphate groups across the base pairing (r = 16 ~ 18 A) is almost negligi- 
ble. The Debye-Huckel interactions is most appropriate for monovalent counterions like Na + . 

Models for linker or "handles": In laser optical tweezer (LOT) experiments the RNA 
molecules are attached to the polystyrene beads by an RNA/DNA hybrid handles or linker. 
A schematic illustration of the pulling simulations that mimic the experimental setups for the 
laser optical tweezer (LOT) (Fig|21-A) is shown in Fig|2}B. Globally, the principles involved 
in atomic force microscope (AFM) and LOT for mechanical unfolding of biomolecules are 
essentially the same except for the difference in the effective spring constant. The spring 
constant of the nearly harmonic potential of the optical trap is in range k = 0.01 — 0.1 
pN/nm whereas the cantilever spring in AFM experiments (FigI2-B) is much stiffer and 
varies from 1 — 10 pN/nm. To fully characterize the RNA energy landscape it is necessary 
to explore a wide range of loading rates 22|| . To vary rj we have used different values of k 



in the simulations. 

We simulate mechanical unfolding of RNA hairpins either by applying a constant force on 
the 3'- end with the 5'- end being fixed (unfolding at constant force), or by pulling the 3'- end 
through a combination of linkers and harmonic spring (FiglljB) at a constant speed in one 
direction (unfolding at constant loading rate). Comparison of the results allows us to test the 
role of the linker dynamics on the experimental outcome. If the linker is sufficiently stiff then 
it should not affect the dynamics of RNA. On the other hand, unfolding at constant loading 
rate (r/ = k X v, where v is the pulling speed and the stretching force (fs) is computed 
using / = k x 5z with Sz being the displacement of the spring) can be modulated either by 
changing k or v. Instead of k an effective spring constant /c e //, with k~L = k^+k^^+k^, 
should be used to compute the loading rate. Typically, k^ ker and k~^ ol (where k mo i is the 
stiffness associated with the model) are small and, hence k e ff ~ k. In this setup, the relevant 
variables are k, v, and the contour length (L) of the linker and the effective stiffness of the 
linker. We explore the effects of these factors by probing changes in the force extension 
curve (FEC). The variations in L are explored only using constant loading rate simulations. 
Manosas and Ritort 0] used an approximate method to model linker dynamics. 

The energy function for the linker molecule is 

jV-l , N-2 

Vl = ^2 y^'^ 1 ~ b "> 2 ~ S • h+l,i+2 (9) 

i=l i=l 

where r iti+ i and fii+i are the distance and unit vector connecting i and i + 1 residue, 
respectively. For the bond potential we set ks = 20 kcal/{mol ■ A 2 ) and b = 5 A. This 
form of the energy function describes the worm-like chain (WLC) [23| (appropriate for 
RNA/DNA hybrid handles used by Liphardt et. al. 0) when Ua is large. The linker 



7 



becomes more flexible as the parameter describing the bending potential, fc^, is reduced. 
Thus, by varying the changes in the entropic elasticity of the linker on RNA hairpin 
dynamics can be examined. We use = 80 kcal/mol or 20 kcal/mol to change the 
flexibility of the linker. For the purpose of computational efficiency, we did not include 
excluded volume interactions between linkers or between the linker and RNA. When 
linkers are under tension the chains do not cross unless thermal fluctuations are larger 
than the energies associated with force. To study the linker effect on force extension 
curves (FEC), we attach two linker polymers each with the contour length L/2 to the ends 
of the hairpin and stretch the molecule using the single pulling speed 0.86 x 10 2 fim/s 
with spring constant k = 0.7 pN/nm. The total linker length is varied from L = (10—50) nm. 

Simulations: We assume that the dynamics of the molecules (RNA hairpins and the 
linkers) can be described by the La ngey in equation. The system of Langevin equations 
is integrated as described before 0, E3 |. Using typical values for the mass of a bead in 
a nucleotide (B{, Si or P,,), m = 100 g/mol ~ 160 g/mol, the average distance between 
the adjacent beads a = 4.6 A, the energy scale e/j = 1 ~ 2 kcal/mol, the natural time is 
t~l — (^^sl) 1 / 2 = 1.6 ~ 2.8 ps. We use tl = 2.0 ps to convert the simulation times into 
real times. To estimate the time scale for thermal and mechanical unfolding dynamics we 
use a Brownian dynamics algorithm [i^, |2^| for which the natural time for the overdamped 
motion is th = ^tl. We use £ = SOr^ 1 , which approximately corresponds to friction 
constant in water, in the overdamped limit. To probe the thermodynamics and kinetics of 
folding we used a number of physical quantities (end-to-end distance (R), fraction of native 
contacts (Q), the structural overlap function (x), number of hydrogen bonds n bond , etc) to 
monitor the structural change in the hairpin. 

Computation of free energy profiles: We adopted the multiple histogram technique 
|27L l28l | to compute the thermodynamic averages of all the observables at various values of 
T and /. For example, the thermodynamic average of the fraction of native contact, Q, can 
be obtained at arbitrary values of T and / if the conformational states are well sampled over 
a range of T and / values. The thermodynamic average of Q is given by 

V n r -<E-fR)/T ELiMg.-R.Q) 

Z^lE,R,Q V e yTR (F k -(E-f k R))/T k 

Q(T, f) = Ek y K k h(ERO) = E W( T > / ( 10 ) 

r -(E-fR)/T l^ k=1 h k {E,R,Q) ' 

l^E,R,Q e Y: k=1 n k e( F k-( B ~fkR))/T k Q 

where K is the number of histograms, hk(E,R,Q) is the number of states between E and 
E + 5E, R and R + 5R, Q and Q + 5Q in the fc-th histogram, rik = J2e r q hk(E, R, Q), 
Tk and fk are the temperature and the force in the simulations used to generate the k— th 
histogram, respectively. The free energy, Fk, that is calculated self-consistently, satisfies 

e -F r /T r = c -(E-f r R)/T r Y.k=l h k{ E i R iQ) / n x 

^ n k e^-i E -fhR))/T h ' 1 ; 

E,R,Q Z^fc=l "''c c 

Using the low friction Langevin dynamics, we sampled the conformational states in the (T,f) 
in the range {0 K < T < 500 K, f = 0.0 pN} and {0.0 pN < f < 20.0 pN, T = 305 K}. 
Exhaustive samplings around the transition regions at {305 K < T < 356 K, f = 0.0 pN} 
and {5.0 pN < /' < 7.0 pN, T = 305 K} is required to obtain reliable estimates of the 
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thermodynamic quantities. The free energy profile, with Q as an order parameter, is given 
by 

F[Q(T, /)] = F Q (T, f) - k B T logP[Q(T, /)]. (12) 
where F D (TJ) = -k B TlogZ(TJ), Z(TJ) = Ee,r,q K E h^ R }% )/Th 

' fc=l 

and P[Q(T, /)] is defined in Eq.fjlO)). The free energy profile F(R) with R as a reaction 
coordinate can be obtained using a similar expression. 



RESULTS 



Mechanisms of thermal denaturation and forced-unfolding are different: We 

had previously reported 0] the thermodynamic characteristics of the P5GA hairpin as a 
function of T and /. The native structure of the hairpin, that was determined using a 
combination of multiple slow cooling, simulated annealing, and steepest descent quench, 
yielded conformations whose average root mean square deviation (RMSD) with respect to 
the PDB structure is about 0.1 A hj. The use of Go-like model leads to a small value of 
RMSD. From the equilibrium (T,f) phase diagram in (TiJ the melting temperature T m « 341 
K at fs = 0. Just as in thermal denaturation at T m at zero / the hairpin unfolds by a weak 
first order transition at an equilibrium critical force f c . Above f c , which is a temperature- 
dependent, the folded state is unstable. 

To monitor the pathways explored in the thermal denaturation we initially equilibrated the 
conformations at T = 100 K at which the hairpin is stable. Thermal unfolding (/ = 0) was 
initiated by a temperature jump to T = 346 K > T m . Similarly, forced- unfolding is induced 
by applying a constant force fs = 42 pN to thermally equilibrated initial conformation at 
T = 254 K Q. The value of f s = 42 pN far exceeds f c = 15 pN at T = 254 K. Upon 



thermal denaturation, the nine bonds fluctuate (in time) stochastically in a manner that is 
independent of one another until the hairpin melts (Fig|3]A). Forced-unfolding, on the other 
hand, occurs in a directed manner. Mechanical unfolding occurs by sequential unzipping 
with force unfolding the bond, from the ends of the hairpin (beginning at bond 1) to the 
loop (FigHB). 

The differences in the folding pathways are also mirrored in the free energy profiles. 
Assuming that Q is an adequate reaction coordinate for thermal unfolding j2^] we find, by 
comparing F(Q) at T = 100 K and T = 346 K, that thermal melting occurs by crossing a 
barrier. The native basin of attraction (NBA) at Q = 0.9 at T = 100 K is unstable at the 
higher temperature and the new equilibrium at Q ~ 0.2 is reached in an apparent two state 
manner (Fig|3]-C). Upon "directed" mechanical unfolding the free energy profile F(R) is 
tilted from the NBA at R ~ 1.5 nm to R ~ 12 nm at which the stretched states are favored 
at / = 42 pN. The forced-unfolding transition also occurs abruptly once the activation 
barrier at R ~ 1.5 nm (close to the folded state) is crossed (FiglUD). 



Free energy profiles and transition state (TS) movements: Based on the prox- 
imity of the average transition state location, Ax^ s , it has been suggested that folded states 
of RNA and proteins 0] are brittle. If the experiments are performed by stretching at a 
constant loading rate then Ax^ s is calculated using /* ~ k B T / Ax^ s log Tf [3(| where /* is 
the most probable unfolding force and r/, the loading rate, is rj = df/dt = kv. Substantial 
curvatures in the dependence of /* on logr^ ([/*,logr^] plot) have been observed especially 
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if Tf is varied over a wide range 22l . Similarly, in constant force unfolding experiments 



Axp 5 is obtained from the Bell equation |15[ that relates the unfolding rate to the applied 
force, \ogku = log k^ + f Ax T ^ jk^T where ky is the unfolding rate in the absence of 
force. In the presence of curvature in the [f*,logrf] plots or when the Bell relationship is 
violated [3lJ it is difficult to extract meaningful values of k^ or Ax^ s by a simple linear 
extrapolation. By carefully examining the origin of curvature in the [f *, log rf] plot or in 
k\j as we show that in the unfolding of hairpin the observed non-linearities are due to 
movements in the transition state ensemble i.e., Ax^ s depends on fs and rj. 

Unfolding at constant loading rate: We performed forced unfolding simulations by 
varying both the pulling speed v and the spring constant k so that a broad range of loading 
rates can be covered. The unfolding forces at which all the hydrogen bonds are ruptured are 
broadly distributed with the average and the dispersion that increase with growing loading 
rates (FigEJA). The plot of /* as a function of logrj (FigEJ-B) shows marked departure 
from linearity. The slope of the plot ([/*,logr/]) increases sharply as r/ increases. There are 
two possible reasons for the increasing tangent. One is the reduction of Ax^ s , which would 
lead to an increase in the slope (/c^T/Ax^ 5 ) of /* vs log 77. The other is the increase of 
curvature at the transition state region, i.e., barrier top of free energy landscape. Regardless 
of the precise reason, it is clear that the standard way of estimating Ax^ s using /* at large 
loading rates results in a very small value of Ax^ s . We estimate Axp from [/*,logr/] plot 
to be 4 A at rj « 10 5 pN/s. The estimated value of Ax^ s is unphysical because 4 A is less 
than the average distance between neighboring P atoms. The minimum pulling speed used 
in our simulations is nearly five orders of magnitude greater than in experiments. The use 
of large loading rate results in small values of Axp S . If the simulations can be performed 
at small values of ry we expect the slope of /* vs logry to decrease, which would then give 
rise to physically reasonable values of Axp at low r/. Our simulations suggest that the 
curvature in the plot of /* as a function of logrj is due to the dependence of Ax^ s on rf 
and not due to the presence of multiple transition states (HI)- As a result, extrapolation to 
low rf values can give erroneous results (Fig0J-B). 

Unfolding at constant force: To monitor the transition state movements we performed 
a number of unfolding simulations by applying fs > 20 pN at T = 290 K. The unfolding 
rates are too slow at fs < 20 pN to be simulated. Nevertheless, the simulations give strong 
evidence for force-dependent movement of Ax^ 3 . For a number of values of fs in the range 
20 pN < fs < 150 pN we computed the distribution of first passage times for unfolding. 
The first passage time for the 2-th molecule is reached if R becomes R = 5 nm for the first 
time. From the distribution of first passage times (for about (50-100) molecules at each fs) 
we calculated the mean unfolding time. Just as for unfolding at constant rf the dependence 
of log tjj on fs shows curvature (FigJSjA), and hence deviates from the often used Bell model 
[HI. By fitting t\j to the Bell formula (tu = Tue~f sAx F S / kBT ), over a narrow range of fs we 
obtain Ax^ s ~ 4A which is too small to be physically meaningful at fs = 0. 

Insights into the shift of Ax^ s as fs increases can be gleaned from the equilibrium 
force-dependent free energy F(R) as a function of R. The one-dimensional free energy 
profiles F(R) show significant movements in Ax^ 3 as fs changes (FigJSJ-C). As fs increases 
AxJ, s decreases sharply, which implies that the unfolding TS is close to the folded state. 
At smaller values of fs the TS moves away from the native state. At the midpoint of 
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the transition Axp s ~ 5.5 ran which is about half-way to the native state. The result 
for Ax^/Ru « 1/2 (i^ is the average value of R in the unfolded state) is in accord 
with experiments |5j which were done at forces that are not too far from the equilibrium 
unfolding force. Although Axp^ is dependent on the RNA sequence it is likely that, for 
simple hairpins Ax F /Ru ~ 1/2. The prediction that Axp is dependent on fs is amenable 
to experimental test. 

Force-quench refolding times depend (approximately) exponentially on fq: 

One of the great advantages of force quench refolding experiments [H is that the ensemble 
of conformations with a predetermined value of R can be prepared by suitably adjusting 
the value of fs 0]. Because force-quench refolding can be initiated from conformations with 
arbitrary values of R, regions of the energy landscape that are completely inaccessible in 
conventional experiments can be probed. In order to initiate refolding by force quench we 
generated extended conformations with R = 13.5 nm, using fs = 90 pN at T = 290 K. 
Subsequently, we reduced the force to fq in the range 0.5 pN < fq < 4 pN. For these values 
of fq the hairpin conformation is preferentially populated at equilibrium. The probability 
Pfy(t) that the RNA hairpin remains unfolded at six fq values decays non-exponentially 
(FigEJ). The mean refolding times 7>(/q), upon force quench, that are computed using 
P u {t) (rp = J™tPu(t)dt), show that in the range 0.5 pN < f Q < 4 pN (FigEJ-B) 

tfUq) = r F exp {fqAx^ 8 /k B T) (13) 

where Axjf is the distance from the unfolded basin of attraction to the refolding transition 
state, and r F is the refolding time in the absence of force. Linear regression gives t f rs 290 
fis and AxJj S ~ 1 nm. The value of Axjf , which is obtained from kinetic simulations, is in 
accord with the location of AxJ obtained directly from the equilibrium free energy profiles 
F(R) (FigGJC). In the 0.5 pN < f Q < 4 pN range the distance from UBA to the TS is 
approximately 1 nm and is a constant. Thus, kinetic and thermodynamic data lead to a 
consistent picture of force-quench refolding. 

If the simulations are done with fq = pN then we find that Tp(fq = 0) ~ 191 lis 
(FigEJ-B) which differs from r F obtained using Eq. (jl3j) . When Jq is set to zero, the 3' 
end fluctuates whereas when /q / the 3' end remains fixed. The rate limiting step in 
the hairpin formation is the trans— ►gauche transitions in the dihedral angles in the GAAA 
tetraloop region (see below). When one end is free to fluctuate, as is the case when fq = 0, 
the trans— >gauche occurs more rapidly than fq ^ 0. The difference between Tp(fq = 0) and 
t f is, perhaps, related to the restriction in the conformational search among the compact 
structures which occurs when one end of the chain is fixed. 

Metastable intermediates lead to a lag phase in the refolding kinetics: The 

presence of long lived conformations (see below), with several incorrect dihedral angles in 
the GAAA tetraloop, is also reflected in the refolding kinetics as monitored by Pu{t). If 
there are parallel routes to the folded state then Pjj(t) can be described using a sum of 
exponentials. The lag kinetics, which is more pronounced as fq increases (see especially 
fq = 4 pN in FigJUJ) can be rationalized using a kinetic scheme S / F where S is 
the initial stretched state, / is the intermediate state, and F is the folded hairpin. Setting 
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P v (t) = P s (t) + = 1 - P F (t), we obtain 

Puit) = (r 2 e-^ 2 - ne^ 1 ) . (14) 

T2 - 7i 

From FigJBJ which shows the fits of the simulated Pjj{t) to Eq.fjl4jl. we obtain the param- 
eters (ti, t 2 ) at each /q (see caption to FigJU] for the values). If folding is initiated by 
temperature-quench t\ <C t 2 so that Pu(t) ~ e - '/ 7 " 2 . Explicit simulations show that thermal 
refolding occurs in a two-state manner (data not shown). In force-quench refolding both 
7~i and r 2 increase as /q increases and ti/t 2 ~ 0(1) at the higher values of fq. There 
are considerable variations in the structures of the metastable intermediate depending 
on fq. The variations in the conformations are due to the differences in the number 
of incorrect or non-native dihedral angles. As a consequence there are multiple steps in 
force quench refolding in contrast to forced-unfolding which occurs in an all-or-none manner. 



Trans^gauche transitions in the GAAA tetraloop dihedral angles lead to long 
refolding times: It is of interest to compare the refolding times obtained from stretched 
ensemble (tf( fq)) and the refolding time (t>(T)) from thermally denatured ensemble. In a 
previous paper [lj], we showed that Tp(fq = 0) = 15t>(T) (Fig|5]-B). The large difference in 
refolding times may be because the initial conditions from which folding commences are vastly 
different upon force and temperature quench 0, 3^]. The fully stretched conformations, 
with R = 13.5 nm, are very unlikely to occur in an equilibrated ensemble even at elevated 
temperatures. The canonical distribution of thermally denatured conformations even at 
T = 1500 K Tp) shows that the probability of populating conformations with R = 13.5 
nm (Fig0-A) is practically zero. Thus, folding from thermally denatured ensemble starts 
from relatively compact conformations. In contrast, the initial condition for force-quench 
refolding can begin (as in our simulations) from fully stretched conformations upon force- 
quench. Both R and the radius of gyration (R g ) undergo substantial changes en route 
to the NBA. Indeed, the refolding trajectories from extended conformations exhibit broad 
fluctuations in R(t) in the order of (25-75A) for long time periods, followed by a cooperative 
reduction in R at the final stage (Fig0-B). 

The long refolding times upon force-quench starting from fully stretched conformations 
may be generic for folding of globular proteins as well. Recent experiments on force-quench 
refolding of poly-ubiquitin (Ub) [H show that R(t) for proteins exhibit behavior similar to 
that shown in Figd-B. The resulting /g-dependent refolding times for poly-Ub (0.1 — 10 sec) 
are considerably larger compared to r F (^i 5 msec ^^,^5^) in the absence of force. Because 
the collapse of a single ubiquitin molecule in solution occurs in less than a millisecond, the 
origin of the long refolding times has drawn considerable attention (3^, HH . The microscopic 
model considered here for RNA can be used to shed light on the origin of the generic long 
refolding times upon force quench. 

From the phase diagram (see Fig. (2) in 01) it is clear that the routes navigated by 
RNA hairpin upon thermal and force quench have to be distinct. While the distinct initial 
conditions do not affect the native state stability (as long as the final values of T and fq are 
the same) they can profoundly alter the rates and pathways of folding. The major reason 
for the long force-quench refolding times in RNA hairpins is that in the initial stretched 
state there is a severe distortion (compared to its value in the native state) in one of the 
dihedral angles. The 19-th dihedral angle (found in the GAAA loop region) along the sugar- 
phosphate backbone is in g + conformation in the native state while in the initial stretched 
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conformation it is in the t state (FigJSJ-A). Thus, if all the molecules are in the fully stretched 
conformation then the 19 th dihedral angle in each of them has to, during the force-quench 
refolding process, undergo the t — > g + transition in the GAAA tetraloop region in order to 
fold. The enthalpic barrier associated with the t — > g + transition is about lksT (FigJHJ-B). 
However, this transition is coupled to the dynamics in the other degrees of freedom and such 
a cooperative event (see Fig|7|-B for examples of trajectories) makes the effective free energy 
barrier even higher. Although significant fluctuations are found in thermally denatured 
ensemble at T = 500 K, they are not large enough to produce non-native dihedral angles 
in the GAAA tetraloop. The dihedral angles in thermally denatured conformations do 
not deviate significantly from their values in the native conformation (FigJSJ-C). In contrast 
upon fully stretching P5GA dihedral angles in the GAAA tetraloop adopt non-native values 
(FigED). 

The time scale for the t — > g + transitions can be inferred from the individual trajectories. 
Typically, there are large fluctuations due to g + <-» t <-» g~ transitions in the dihedral 
angles in the flexible loop region (i — 19 — 24). For the trajectory in FigUB we observe 
the persistence of incorrect dihedral angle in the loop region for t ~ 300 fis. At t > 300 
lis the native-like dihedral angles form. Subsequently, the formation and propagation of 
interaction stabilizing the native RNA hairpin takes place. These dynamical transitions 
are clearly observed in Figs El- B and |UJ-C. The observed mechanism is reminiscent of a 
nucleation process. We conclude that the formation of the flexible loop with all the dihedral 
angles achieving near native values is the rate limiting step in the refolding kinetics of RNA 
hairpins upon force-quench starting from fully stretched state. It should be stressed that 
the rate limiting step for thermal refolding of P5GA or force-quench refolding starting from 
partially stretched conformations (see below) is different. The observation that the zipping 
of the hairpin takes place upon synchronous formation of all the native-like dihedral angles 
suggests the presence of a high entropic barrier. The crossing of the entropic barrier results 
in slow refolding if P5GA is fully stretched. 

Linker effects on RNA force-extension curves : In LOT experiments the force- 
extension curves (FECs) are measured for the handle (H)-RNA-handle(H) construct E[. 
To unambiguously extract the FEC for RNA alone (from the measured FEC for H-RNA-H 
construct) the properties of the handle, namely, the contour length Lh and the persistence 
length ljf cannot be chosen arbitrarily. In the experiments by Liphardt et. al. Lh = 320 nm 
and If = 50 nm for the RNA/DNA hybrid handle. We expect that both A = If /lf NA and 
lp will affect the FEC curves of the object of interest, namely, the RNA molecule. To discern 
the signature for the force-induced transition in the RNA hairpin alone from the FEC for 
H-RNA-H construct A has to be large. If we assume lf NA ~ 1 nm then the experimental 
value of A ~ 50 which is large enough to extract the transitions in the RNA hairpin. If A ~ 1 
then the entropic fluctuations in the handle can mask the signals in RNA j3(| . Similarly the 
end-to-end distance fluctuation in the handle, 5R, should be smaller then the extension in 
RNA. Because 5R grows with Lh (see below) it follows that if very long Lh is used even 
with A ^> 1 the signal from RNA can be masked. The square of the fluctuations in the 
end-to-end distance 5R of the linker is given by {SR) 2 = d(R) /d(j3fs) where (R) is the mean 
end-to-end distance. For WLC, that describes the linkers, we expect (R) ~ {2L H lp) l l 2 for 
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small fs and (R) ~ L H for large fs provided L H is long. Thus, 



5R 



(L H l«y/\k B T/f s )^ (f s < k B T) 
Lf(k B T/f s )^ (f s » k B T) 



(15) 



The mean fluctuation in the extension of the spring is, using equipartition theorem, given 



In order for the signal from RNA to be easily discerned from the experimentally measured 
FEC, the expansion of end-to-end distance of the molecule at transition should be larger than 
SR and 8x . Since 5R grows sub linearly with the linker length (Eq.([15j)) the attachment of 
large linker polymer can mask the transition signal. These arguments show the characteristics 
of the linker can obscure the signals from RNA. 

The FECs in the simulations can also be affected by non-equilibrium effects due to the 
linker dynamics. In the handle(H)-RNA-handle(H) construct considered here the force ex- 
erted at one of the linker depends on the angle between fs and the end of the linker. The 
initial event in force transmission along the contour of the H-RNA-H construct is the align- 
ment of the molecule along the force direction. The characteristic time for force to reach 
RNA so that unfolding can occur is f c /vf where f c is the critical force. Non-equilibrium ef- 
fects due to linker dynamics become relevant if tr, the time scale for alignment of H-RNA-H 
along the force direction, tr > f c / r f- This condition is not relevant in experiments which are 
conducted at small values of r/. However, it is important to consider non-equilibrium effects 
in simulations which are performed at high ry values. The characteristic time depends on 
L H and 

We validate these arguments by obtaining FEC for H-P5GA-H by varying L B and A = 
lp /lp NA - Using the WLC for the linker we calculated FEC where the L B is varied from 
(10 — 50) nm. In order to observe rapid unfolding we have carried out our simulation at 
the pulling speed v = 0.86 x 10 2 fim/sec and k = 0.7 pN/nm. Under these conditions 
non-equilibrium effects are relevant for the linker [23| which is not the case in experiments. 
The FECs show clearly a plateau in the range 20 pN < fs < 40 pN which corresponds to 
the two-state hairpin opening (FigllU|). For the experimentally relevant plot (FigllOt-A) that 
shows FEC for H-P5GA-H the transition plateau is present at all values of L B . However, 
when L B reaches 50 nm the signal from P5GA is masked. The FEC for P5GA alone fFigfTUl 
B) shows modest increase in the unfolding force as L B increases. Similarly, we find the value 
of unfolding force also increases as the linker flexibility increases. These observations are 
due to non-equilibrium effects on the linker dynamics because of the relatively large values 
of rj used in the simulations. 

Our simulations show that, at high loading rates, the length of the handle is also 
important. This issue is not relevant in experiments in which loading rates are much 
smaller. However, they become important in interpreting simulation results. In our case 
rf(= kv) is 6 x 10 4 times that used in experiments. At such high rf non-equilibrium effects 
control the linker dynamics [23j]. Thus, to extract unfolding signatures from RNA alone it 
is necessary to use high value of A and relatively short values of L. In other words FEC, 
when rj is varied, may be a complicated function of lp /lp NA and L B /lp ■ 



by 




(16) 
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Force-quench refolding of P5GA with attached linkers: It is convenient to monitor 
force-quench refolding of RNA alone using simulations. A similar experiment can only be 
performed by attaching handles to RNA. In such an experiment, which has not yet been 
done (however, see Note added), the H-RNA-H would be stretched by a stretching force 
fs > fc so that RNA unfolds. By a feedback mechanism the force is quenched to /q ^ 0. 
We have simulated this situation for the H-P5GA-H construct with Lh = 15 nm, lp = 30 
nm for each handle. We chose very stiff handles (Lh < lp) so that the dynamics of RNA 
can be easily monitored. The end-to-end distance of the H-P5GA-H system as a function 
of t with fs = 90 pN and /q = 2 pN shows a rapid decrease from R sys = 44 nm to about 
Rsys = 37 nm in less than about 100 fis (FigJTTJ). The use of large values of fs will not affect 
the results qualitatively. The value of R sys fluctuates around 36 nm for a prolonged period 
and eventually R sys attains its equilibrium value around 30 nm. 

Upon decomposing R sys into contributions from P5GA and the handle we find that the 
major changes in R sys occur when RNA undergoes the folding transition (compare the top 
and middle panels in Fig fTT|) . The time dependence of Rh, which monitors only the dynamics 
of the linker, shows that with /q = 2 pN after the initial relaxation Rh fluctuates around 
its equilibrium value (bottom panel in Fig [TTj) . 

From these simulations and others for different /q values we can make a few general 
comments that are relevant for experiments, (a) As long as A(= lp /lp NA ) is large enough 
the qualitative aspect of RNA folding can be obtained from the dynamics of R sys alone. 
However, if Lh/1? S> 1 then large transverse fluctuations of the linker can interfere with the 
signal from RNA molecule, (b) To obtain quantitative results for the dynamics of RNA (i.e., 
Rm as a function of t) the dynamics of the handle upon force relaxation has to be described 
accurately. Upon fs — > /q quench the dynamics of the handle cannot be described using 
Langevin equation using the equilibrium force. Instead, the relaxation behavior must be 
determined by solving the Langevin equation for the WLC energy function j3?J which is 
subject to fs — > /q quench. 



DISCUSSIONS 



Transition state movement and Hammond postulate for force: Hammond 
postulate [sH 3^] is widely used to qualitatively predict the nature of transition state in 



the chemical reactions of organic molecules. In recent years a number of protein folding 
experiments have been interpreted using generalization of the Hammond postulate . The 
Hammond postulate states that if a transition state and an unstable intermediate, occur 
consecutively during a reaction process and have nearly the same energy, their interconversion 
will involve only a small reorganization of the molecular structure (3^ |. In the context of RNA 
folding, the Hammond postulate suggests that the position of the transition state along 
the reaction coordinate is shifted towards the destabilized state, either folded or unfolded, 
depending on the nature of perturbation. The Hammond behavior is most vividly seen in the 
free energy profiles F(R) (FigJSJ-C) that pictorially describe mechanical unfolding of RNA 
hairpins. As / is increased the unfolded state is preferentially stabilized. From the Hammond 
postulate we would infer that the major TS should be more native-like as / increases. The 
force-dependent F(R) as a function of R indeed confirms (FigJSJ-C) that as / increases Ax^ s 
becomes closer to the native folded hairpin conformation. 

RNA hairpins also denature upon heating. To ascertain the variation in the location of 
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the TS as temperature is changed we have calculated the free energy F(Q) as a function 
of Q at several values of T at / = (FiglT^I-A). Although the location of the TS follows 
Hammond behavior there is very little change in the TS ensemble over the temperature range 
examined. Thus, the changes in the TS are very dramatic when unfolding is induced by force 
compared to thermal denaturation. 

Hammond behavior can be quantified using the Leffler's proportionality constant ct x which 
measures the energetic sensitivity of the transition state relative to the native states when the 
population shift is induced by a perturbation x (111 EH- For mechanical unfolding (x = fs) 

dAF*(R)/df s Ax T F s 
af dAF UF (R)/df s Ax UF 1 } 

Using the free energy profile we computed «j as a function of AF UF or / (FigIT2T-B). The 
shift in the transition state is, quantified by ay in the range < af < 1, and the shift 
rate (or self- interaction parameter, pj = dctf /OAFuf jlH) has maximum in the force range 
4 < fs < 10 pN. As AFyp decreases (the UBA is stabilized with respect to the NBA) af 
decreases, which implies that the TS becomes increasingly native-like fFigfO-B). The inset 
in FigfT2l-B shows dramatically the changes in af with respect to fs- The largest changes in 
af occurs as fs approaches the T-dependent (T = 290 K) f c ~7 pN. A similar plot of a? 
as a function of T shows practically no change in ax- From this analysis we conclude that 
the nature of the transition state ensemble is different in mechanical unfolding and thermal 
denaturation. 

The transition state movement with force is very sensitive to the shape of the barrier in 
the vicinity of the transition state. The free energy profile near the barrier (x ~ x ts ) can be 
expanded as F(x) ~ F(x ts ) — ^F"(x ts )(x—x ts ) 2 + - ■ ■■ Upon application of the stretching force 
F(x) is tilted by —fs ■ x. The new barrier position (xf s EW ) is at xf s EW ~ x ts — fs/F"(x ts ). 
For a sharp transition barrier (x ts F"(x ts ) 3> fs), the force will not affect the position of 
the transition state {x^ s EW ~ Xt s )- If the transition barrier is broadly distributed as in the 
unzipping pathway of RNA hairpins, the structure of transition state progressively changes 
as the magnitude of the force is varied. Typically, folding transition states are shallow and 
broad. As a result in biomolecular folding or unbinding, which involve formation or rupture 
of non-covalent interactions, we predict that the location of the TS depends on fs and 
temperature. The assumption of a fixed TS used to interpret 0, H2, H3| experimental results 



is not valid. In addition, sequential |22[ and/or parallel pathways to the stretched state 
transition [3, EH are also possible. These observations suggest that a careful inspection 
of nonlinearity in the Arrhenius plot and af will be required to unravel barriers to unfolding. 

Entropic barriers and long refolding times from fully stretched state: We 

have proposed that the long refolding time in P5GA upon force-quench from the initial 
stretched conformations is due to entropic barriers. The rate limiting step in the force-quench 
refolding of P5GA is the t — > g + transition in the nucleotides near the GAAA tetraloops. 
We analyze the simulation results by adopting a model proposed by Zwanzig [46[. In this 
Ising-like model each degree of freedom (in our case the dihedral angle) is presumed to exist 
in the "correct" state (native) and "incorrect" or non-native state. Suppose that the energy 
difference between the incorrect and correct states is e and that there are correct dihedral 
angles required for forming the hairpin loop. Let the free energy of loop stabilization be ei oop . 
The energy, E n , of a conformation with n-incorrect dihedral angles is E n = ne — S n o€i oop . If 
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we assume that the dihedral angle is a discrete variable with l + v states [y is the number 
of incorrect states) then the partition function is 



JV 

E 

n=0 



Z N = V ( N ) ^ e -P(neS n0 ei oop ) = e (3e loop + (1 _|_ ve ~^)N _ \ 



n 



The thermal probability of realizing a conformation with n incorrect dihedral angles is 

N \ u n e -(J(ne-8 n0 e loop ) 



Pn(eq) = ^ . (19) 

Zjn 

The free energy profile, with n playing the role of a reaction coordinate for dihedral angle 
transitions, is AF(n) = ne - 5 n0 e loop - k B T \ogv n ( N n ) . For P5GA N = 6, e « 2k B T, and 
e loop « 7.6£; B T (= V$%gf£ 4Bl3 + V^ Bli = 5.1k B T + 2.5k B T). The free energy profile AF(ra) 
(Fig 113]) shows that the barrier depends only weakly on v. Because the dihedral angle is a 
continuous variable, we use v as an undetermined parameter. For v = 13 the free energy 
barrier (AF F ) is ~ 2.7k B T, which leads to the observed increase (tf(/q = 0) = 15t>(T)) in 

the refolding time by factor of 15 (= e AF f/ fcsT ). 

The entropic barrier ~ (5 — 6)k B T is significantly larger than the free energy barrier AFfj 
in the absence of force. The barrier to the formation of conformations with (tttttt) state in 
the loop region is large enough that they do not form by thermal fluctuations, and hence are 
irrelevant when refolding is initiated by temperature quench. However, such conformations 
are populated with near unit probability when fully stretched by mechanical force. When 
folding is initiated by force-quench from extended conformations (like conformation I in 
FigJSJ-E), metastable conformations with incorrect dihedral angles in the loop regions are 
formed. These are characterized by plateaus in the dynamics of R (FigJUJ-A). The crossing 
of the entropic barrier that places the loop dihedral angles in the native-like gauche state 
results in the slow refolding of P5GA hairpins. Once the loop is formed the zipping process 
quickly stabilizes the hairpin so that barrier crossing in reverse direction is unlikely to occur 
at low forces. 

To further validate the proposed mechanism we performed simulations in which the value 
of the initial stretching force fs is not large enough to fully extend the P5GA hairpin. In these 
simulations, the ensemble of initial structures is prepared so that they contain the preformed 
GAAA loop with only single bond before the loop region that is intact (see conformation II in 
Fig|SJ-E). The refolding kinetics follows a single exponential decay with a mean refolding time 
of ~ 33 fis, which is only 2-3 times longer than the refolding time of thermally denatured 
states. This value is much shorter than the refolding time from the fully stretched states 
(> 191 [is, see FigUB). These simulations also show that rp is a function of both fs and 
/q- 

A byproduct of this analysis is that the appropriate reaction coordinate in force-quench 
refolding of RNA hairpins may be a local variable. In the formation of P5GA hairpin the 
local dihedral angles are the relevant reaction coordinates. The local dihedral coordinates, 
that describe the rate-limiting steps in the UBA— >NBA transition, are hidden in the global 
coordinates such as Q or R. Indeed, there is no correlation between the formation of 
native dihedral angles in the GAAA tetraloop and global order parameters. We infer that 
to describe folding, especially of RNA, multiple reaction coordinates that describe the 
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hierarchical assembly are required. 



Difficulties in extracting energy landscape parameters from single molecule 
force spectroscopy: Several studies 31, 4^] have pointed out the inherent ambiguities in 



quantitatively characterizing the energy landscape from measurable quantities in dynamic 
force spectroscopy. From the plots [/*,logr/] one cannot unambiguousl y o btain the location 
of the transition state(s) or even the number of free energy barriers |42j. The significant 
curvatures in the [/*,logry] plots are usually interpreted in terms of multiple transition states 
(22! 42T |. In our example, the P5GA hairpin unfolds upon application of force by crossing 



a single free energy barrier. Explicit equilibrium F(R) profiles (FigfSJ-C) and experiments 
that have monitored hopping dynamics in P5ab hairpin show that there is only one free 
energy barrier in these simple structures. Nevertheless, [J*,logr/] plot is highly nonlinear 
(FigEJ-B). In the hairpin case we have shown that the non-linearity is due to the dramatic 
changes in Ax^ s as r/ is varied. The standard assumption that AxJj S is a constant breaks 
down, and is likely to be an even more of a severe approximation for RNA with tertiary 
structures. 

To further illustrate the importance of transition state movements we consider a trivial 
one dimensional potential 

E(x) = -eexp(-£x). (20) 

In this barrierless potential a particle is "unbound" if \E(xt s )/kBT\ < 1 where Xt s is the TS 
location. Upon application of a constant f$ the potential becomes 

E{x) = -eexp{-&)-fsx. (21) 
The location of the transition state in the force range < f§ < £e is 

x™ = -hog fs/£e. (22) 

If fs > £e then xJj S (f) = and the particle is always unbound (FiglTH-A). The changes 
in xJj S can lead to significant deviations from the Bell equation even in constant f$- 
experiments. The distribution of unbinding forces upon deforming the potential at constant 
loading rate (fs(t) = Tft) can be analytically obtained (see Eq.(8) in [47[). The [/*,logr/] 
plot from these calculations show non-linearities (FigIT4TB) that are similar to what is 
found for P5GA (FigEJ-B). In both instances the reason for curvature is entirely due to 
rj-dependent changes in the TS location. 



CONCLUSIONS 



We have systematically investigated forced-unfolding and force-quench refolding of RNA 
hairpins. Using a general minimal model for RNA we have obtained a number of new results 
that give a molecular picture of unfolding and refolding of RNA hairpins triggered by force. 
Although they were obtained specifically for P5GA we expect the conclusions to be valid for 
other RNA sequences as well. The specific predictions, that are amenable to experimental 
test, of our study are listed below. 

1) Besides probing the energy landscape using / as a "denaturant" one of the goals of single 
molecule studies is to extract intrinsic parameters like folding (k° F ) and unfolding rates 
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(ky) and the nature of the transition state in the absence of force. However, extraction 
of the kinetic parameters from dynamic force spectroscopy is fraught with difficulties 
because several models can produce similar [/*,logr^] profiles. Here we have shown, for 
a system that has only a single barrier to unfolding, that AxJj S depends dramatically 
on fs and 77. The movements in the TS is intrinsic to properties of RNA hairpins 
and are clearly reflected in the free energy profiles. Combining the present results 



and our previous study [14j we surmise that Ax^ s is dependent on T and /. Thus, 



extrapolations to zero force to obtain reliable estimates of unfolding rates requires not 
only accounting for free energy barriers [3lJ but also on the dependence of AxJj S on T 
and /. Only by performing multiple experiments (or simulations) over a range of fs 
and T can the free energy landscape be fully characterized. 

2) An important prediction of our work is that refolding times upon force-quench 7>(/q) 

from stretched states are much greater than those obtained by temperature quench. 
More generally, 7>(/q) depends, sensitively on the initial value of the stretching force. 
The microscopic origin of the long force-quench refolding times in P5GA has been 
traced to the time needed for the trans— >gauche transition in the GAAA tetraloop 
region. From this observation we predict that refolding time 7~f(/q) should be very 
long compared to thermal refolding times for P5ab which also has the GAAA tetraloop. 
Because the refolding times for RNA hairpins are determined by the local structural 
features in the initial stretched states we suggest that, at a fixed temperature, t>(/q) 
might depend upon only weakly on the helix length or the precise sequence (percentage 
of GC for example). 

3) Dissecting the folding mechanism of RNA is difficult because of an interplay of a number 

of factors JHI . We predict that refolding mechanisms (pathways and the nature of the 
transition state ensemble) by temperature (or by increasing cation concentration) and 
force quench have to be drastically different. In the former, the transition to the low 
entropic NBA proceeds from a high entropy relatively compact state whereas in the 
latter it occurs from a low entropy stretched state [4£|. The predicted dramatic differ- 
ences in the folding mechanisms can be established by probing force-quench refolding 
at fixed T and counterion concentration. 

The present model has a number of limitations. The use of Go model for force-unfolding 
may not be a serious approximation because unfolding pathway are largely determined 
by the native topology |50j]. However, the neglect of non-native interactions will have 
dramatic effect on refolding. At a minimum, the roughness (5e) of the energy landscape 
is underestimated by the Go model. The importance of Se can be assessed by doing 
forced-unfolding experiments over a range of temperature (47l . lo"ll | . Finally, the electrostatic 
interactions in RNA have been modeled in the simplest manner that is only appropriate for 
monovalent cation (H^ j. To address the effect of counterion (Mg 2+ or polyanions) explicit 
modeling of the cations will be required. 

APPENDIX I 

In this appendix we describe the procedure for determining the persistence length of the 
linkers used in our simulations. For the linker molecules, whose energy function is given 
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by Eq.Q, we calculated the persistence length by fitting the worm-like chain end-to-end 
(R) distribution function (Pwr.n(R)) to the simulated P(R)- We adopted Monte Carlo 
simulation with Pivot algorithm [54J to generate a large number of equilibrium conformations 
of the linker molecule. Given ks{= 20 kcal/(mol ■ A 2 )), kA{= 20 kcal/mol or 80 kcal/mol), 
and varying N the number of monomers in the WLC linker, we obtained the unknown 
parameters, namely, the contour length (L) and the persistence length (l p ) by fitting P(R) 
to 

L[l - (R/L) 2 f/ 2 eXP[ ~A{l-{R/L) 

where t = L/l p . The normalization constant C = l/[7r 3 / 2 e~ a a~ 3 / 2 (l + 3a _1 + 15/4a -2 )] 
with a = St/4, satisfies dRPwLc(R) — 1- The dependence of the persistence length of 
the linkers, as a function of N is displayed in FiglTol The quality of the fit improves as N 
becomes larger (data not shown). We also computed the persistence length and the contour 



wlc(R) ~ Th _ r p/ n2 i9/2 ex P[ _ 77i /p/rwJ' ( 23 ) 



nm 



length of P5GA at T > T m using the same fitting procedure, which gives l p ~ 1.5 
and L = 12.5 nm fFiglTBT-inset on the bottom). In ^ur simulations A = l p /lp NA ranges 
from 10 < A < 70. The experimental value of A ~ 50 
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Note added: While the present paper was under review refolding upon force-quench 
of TAR RNA was reported |55[. In accord with the present and our previous studies (lij . 
force-quench refolding times are relatively long. The distributions of refolding times similar 
to the curves in FigEl as a function of /q were not reported in j55[. Thus, it is unclear if 
there is a lag phase in the force quench refolding of TAR RNA. 
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FIGURE CAPTION 



Figure^: Coarse-grained representation of a RNA using three (phosphate (P), sugar 
(S) and base (B)) interaction sites per nucleotide. On the left we present the secondary 
structure of the 22-nt P5GA hairpin in which the bonds formed between base pairs are 
labeled from 1 to 9. The PDB structure [T3| and the lowest energy structure obtained with 
the coarse-grained model are shown on the right. 

Figure[21: A. Schematic illustration of laser optical tweezer (LOT) setup for RNA stretch- 
ing. Single RNA molecule is held between two polystyrene beads via molecular handles with 
one of the polystyrene beads being optically trapped in the laser light. The location of the 
other bead is changed by manipulating it by a micropipette. The extension of the molecule 
through the molecular handles induces the deviation in the position of the polystyrene bead 
held in the force-measuring optical trap. B. Both LOT and AFM can be conceptualized as 
schematically shown. The RNA molecule is sandwiched between the linkers and one end of 
the linker is pulled. The spring constant of the harmonic trap in LOT (or the cantilever in 
AFM experiments) is given by k and v is the pulling speed. 

Figure El: Unfolding pathways upon temperature and force jump. A. The time depen- 
dence of rupture of the bonds is monitored when the temperature is raised from T(= 100 
K) < T m (« 341 K) to T{= 346 K) > T m . The set of nine bonds are disrupted stochastically. 
B. In forced unfolding bonds rip from the ends to the loop regions in an apparent staircase 
pattern. For both A and B the scale indicating the probability of a given bond being intact 
is given below. C. Free energy F(Q) as a function of Q. The stable hairpin with Q « 1 
at T = 100 K becomes unstable upon rapid temperature jump to T = 346 K (blue *). 
Subsequent to the T-jump the hairpin relaxes to the new equilibrium state by crossing a 
small free energy barrier 0.5 kcal/mol) with Q 0.2. The inset shows the equilibrium 
free energy profile at T = 346 K and / = 0. D. Deformation of the free energy profile upon 
application of force. The R dependent free energy F(R) = F(R; fs = fy — fs'R favors the 
stretched state at R — 12 nm when / = 42 pN (see inset). The activation barrier separating 
the UBA and NBA is around (1-2) kcal/mol. 

Figure 0] : Constant loading rate force unfolding. A. The unfolding force distributions 
at different pulling speeds with hard (k = 70 pN/nm, up) and soft springs (k = 0.7 pN/nm, 
down). For the hard spring, the pulling speeds from right to left are v = 8.6 x 10 4 , 8.6 x 10 3 , 
8.6 x 10 2 fim/s. For the soft spring, the pulling speeds are v = 8.6 x 10 4 , 1.7 x 10 4 , 8.6 x 10 3 , 
8.6 x 10 2 , 8.6 x 10 1 fim/s from right to left force peaks. The peak in the distributions which 
are fit to a Gaussian is the most probable force /*. B. The dependence of /* as a function 
of the loading rates, rf. The results from the hard spring and soft spring are combined using 
the loading rate as the relevant variable. The inset illustrates the potential difficulties in 
extrapolating from simulations at large r/ to small values of 77. 

Figure |5] : Kinetics of forced- unfolding and force-quench refolding: A. Plot of force- 
induced unfolding times (tjj) as a function of the stretching force. Over a narrow range 
of force Tu decreases exponentially as / increases. B. Refolding time Tp as a function 
of Jq. The initial value of the stretching force is 90 pN. By fitting Tp using tp(Jq) = 
r£exp (fgAxp/kBT), in the range of 0.5 pN < fq < 4 pN, we obtain Axp w 1 nm and 
t f ~ 290 fis. C. Changes in the equilibrium free energy profiles at T = 290 K F(R) as 
a function of the variable R. We show F(R) at various f$ values. For emphasis, the free 
energies at fs — and at the transition midpoint f$ = 7.5 pN (dashed line) are drawn in 
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thick lines. 

Figure 0: Time dependence of the probability that RNA is unfolded upon force quench. 
In these simulation T = 290 K, and the initial stretching force f$ = 90 pN and Jq (values are 
given in each panel), the quench force, is varied. The simulation results are fit using Eq. fTH) 
which is obtained using the kinetic scheme 5* — — > I F. Here, I represents conformations 
with in certain fraction of incorrect dihedral angles. The time constants (ti,T2), in us, at 
each force are: (81.2, 101.3) at f Q = pN, (159.5, 160.8) at f Q = 0.5 pN, (180.0, 174.8) at 
f Q = 1 P N, (237.6, 240.5) at f Q = 2 pN, (326.8, 335.6) at f Q = 3 pN, and (347.7, 329.7) at 
f Q =}pN. 

Figure [7| : A. The equilibrium distribution of the end-to-end distance at extremely 
high temperature (T = 1500 K). Even at this elevated temperature the fully stretched 
conformations of R = 13.5 nm (arrow) is not found in the ensemble of thermally denatured 
conformations. B. Refolding is initiated by a force quench from the initial value fs = 90 pN 
to /q = 4 pN. The five time traces show great variations in the relaxation to the hairpin 
conformation. However, in all trajectories R decreases in at least three distinct stages that 
are explicitly labeled for the trajectory in green. 

Figure |S] : A. The dihedral angles of the P5GA hairpin in the native state. All the 
dihedral angles are in the trans form except 19-th position of dihedral angle which is in 
the gauche(+) conformation (indicated by orange circle). B. The dihedral angle potentials 
for trans (top) and gauche(+) form (bottom) are plotted using Eq.(j3J). The red lines show 
the potentials in the loop region. C, D. The average deviation of the z-th dihedral angle 
relative to the native state is computed using the 100 different structures generated by high 
temperature (T = 500 K) (C) and by force (R = 13.5 nm) (D). To express the deviations of 
the dihedral angles from their native state values we used 1 — cos (0j — 0°) for z-th dihedral 
angle 0j where 0° is the z-th dihedral angle of the native state. E. A snapshot of fully 
stretched hairpin (I). Note the transition in 19-th dihedral angle undergoes g + — > t transition 
when hairpin is stretched. An example of a partially stretched conformation (II) with the 
GAAA tetraloop and bond 9 intact. Refolding times starting from these conformations are 
expected to be shorter than those that start from fully stretched states. 

Figure |5] : A. A sample refolding trajectory starting from the stretched state. The 
hairpin was initially unfolded to a fully stretched state and /q was set to zero at t ~ 20 /j,s. 
End-to-end distance monitored as a function of time shows that refolding occurs in steps. B. 
The deviation of dihedral angles from their values in native state as a function of time. The 
large deviation of the dihedral angles in loop region can be seen in the red strip. Note that 
this strip disappears around t « 300 us, which coincides with the formation of bonds shown 
in C. /b is the fraction of bonds with pink color indicating that the bond is fully formed. 

Figure El: The force extension curves (FECs) of RNA hairpin at constant pulling speed 
and varying linker lengths and flexibilities. The pulling speed, v = 0.86 x 10 2 /im/s. The 
spring constant, k = 0.7 pN/nm. FECs from different linker lengths are plotted in black (10 
nm), red (20 nm), green (30 nm), blue (40 nm), and orange (50 nm). A. This panel shows 
the experimentally relevant plots, namely, FECs for H-RNA-H construct. The signature for 
the hairpin opening transition region is ambiguous at large values. B. FECs only for 
P5GA corresponding to different Lh values. The FEC for the linker is subtracted from (A). 
The gradual increase of rupture force is observed as L increases. The value of fc^ (Eq.Q) 
of linker polymer used in (A) and (B) is 80 kcal/(mol ■ A). C. Comparison between FECs 
with Lh{= 40 nm) fixed but at different k& values. Red curve is the average of 7 individual 
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FECs that are shown in orange (80 kcalj (mol-A)). Black curve is the average of 9 individual 
FECs in grey (20 kcal/(mol -A)). Less stiff linker leads to slightly larger unfolding force. It 
should be stressed that for both values of k^ the A ratio is large. 

Figure : Refolding trajectory of RNA that is attached to the handles. Linkers with 
Lh = 15 nm and l p = 30 nm are attached to both side of RNA hairpin. The initial force 
(90 pN) stretched P5GA to 14 nm. The value of /q = 2 pN. The top panel shows the 
dynamics of the end-to-end distance, R sys , of H-P5GA-H. The middle panel corresponds 
to the end-to-end distance of P5GA. The end-to-end distances, Rh, of the 3' and 5'-side 
handles fluctuates around its equilibrium value of 14 nm after the initial rapid relaxation. 
The right panels show the decomposition of Rh into the longitudinal and the transverse 
components. The value of R H w 13 nm agrees well with the equilibrium value obtained by 

solving f Q = kgT/lp (r^/Lh + 1/4(1 - f|) 2 - 1/4) with f Q = 2pN. 

FigurelT2l: A. Free energy profiles of Q at different temperatures. Note that the positions 
of the transition states over the temperature variation almost remain constant. B. The 
movement of transition state measured in terms of the Leffler parameter (Eq. (|17Jl ). The 
structural nature of transition state is monitored by the free energy difference between NBA 
and UBA when / is a external variable. The inset shows the variation of a with respect to 
force. The value of T = 290 K. 

Figure EH : The free energy as a function of the number of incorrect dihedral angles 
calculated using the Zwanzig model (see text for details). As the number of distinct values 
(y) that the dihedral angle can take increases the entropic barrier increases. For P5GA 
v = 13 provides the best fit of the model to simulations (see text) . 

Figure El A. Sketch of the one-dimensional potential E(x) as a function of x for several 
values of f$. The transition state location is obtained using E'(xt s ) = 0. The boundary 
separating bound and unbound states is given by \E(x ts )/kBT\ = 1. B. Dependence of the 
most probable "unbinding" force /* as a function of 77. The [/*,logr/] plot for the artificial 
potential is similar to that shown in FigHJ-B. 

Figure : Persistence length l p as a function of contour length L for linkers at 290 K. 
The number of monomers N (Eq.flSJ)) is also shown. The values of l p and L are obtained by 
fitting the end-to-end distance distribution function P(R) generated by simulations to the 
theoretical expression based on a mean field model (Eq.([23|)). An example of such a fit for 
a linker with N = 50 using two different values of k^ (see Eq.fjHJ) is shown in the inset on 
the top. For large N, l p converges to a constant value. The fit of simulated P(R) for P5GA 
computed at T(= 500 K) > T m is shown in the inset on the bottom. From the WLC fit we 
obtain l p « 1.5 nm and L = 12.5 nm for P5GA. 
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